home *** CD-ROM | disk | FTP | other *** search
/ NeXT Education Software Sampler 1992 Fall / NeXT Education Software Sampler 1992 Fall.iso / Mathematics / Notebooks / SigProc2.0 / Packages / SignalProcessing / Digital / InvZTransform.m < prev    next >
Encoding:
Text File  |  1992-08-18  |  38.1 KB  |  1,118 lines
  1. (*  :Title:    Inverse z-Transforms  *)
  2.  
  3. (*  :Authors:    Brian Evans, James McClellan  *)
  4.  
  5. (*
  6.     :Summary:    One-dimensional and separable multidimensional
  7.         inverse z-transforms
  8.  *)
  9.  
  10. (*  :Context:    SignalProcessing`Digital`InvZTransform`  *)
  11.  
  12. (*  :PackageVersion:  2.7    *)
  13.  
  14. (*
  15.     :Copyright:    Copyright 1989-1991 by Brian L. Evans
  16.         Georgia Tech Research Corporation
  17.  
  18.     Permission to use, copy, modify, and distribute this software
  19.     and its documentation for any purpose and without fee is
  20.     hereby granted, provided that the above copyright notice
  21.     appear in all copies and that both that copyright notice and
  22.     this permission notice appear in supporting documentation,
  23.     and that the name of the Georgia Tech Research Corporation,
  24.     Georgia Tech, or Georgia Institute of Technology not be used
  25.     in advertising or publicity pertaining to distribution of the
  26.     software without specific, written prior permission.  Georgia
  27.     Tech makes no representations about the suitability of this
  28.     software for any purpose.  It is provided "as is" without
  29.     express or implied warranty.
  30.  *)
  31.  
  32. (*  :History:    *)
  33.  
  34. (*  :Keywords:    z-transform, region of convergence  *)
  35.  
  36. (*  :Source:    *)
  37.  
  38. (* 
  39.     :Warning:    Keep conditional clauses on the same line; separating them
  40.           by RETURNS will confuse Mathematica (e.g., the package
  41.           context is never reset to Global); you can place a
  42.           newline after an operator (like &&)
  43.  
  44.         The n argument of myinvz must always be a symbol NOT
  45.           an expression (see postinvzrules and local shift operator)
  46.  *)
  47.  
  48. (*  :Mathematica Version:  1.2 or 2.0  *)
  49.  
  50. (*  :Limitation:  *)
  51.  
  52. (*
  53.     :Discussion:  1-D Local state rules base has 62 rules:
  54.             I.   multidimensional hooks     2 rules
  55.             I.   rational transforms       26 rules
  56.             II.  non-rational transforms    9 rules
  57.             III. transform properties      12 rules
  58.             IV.  transform DSP structures   4 rules
  59.             V.   transform strategies       9 rules
  60.  
  61.     Note that  the  partial  fractions  decomposition  strategy is
  62.     implemented  as  two  rules --   one in the rational transform
  63.     pairs section and one in the properties section.
  64.  
  65.     At each step in the inverse rules base, the current expression
  66.     has a local state associated with it.   This state consists of
  67.     a list of nine boolean values.   Each boolean value is associ-
  68.     ated with a strategy.   If  an  element  is  True,  then  that
  69.     strategy has not been tried yet;  if False, then that strategy
  70.     has already been tried, and it will not be tried again.  Thus,
  71.     local state is used  to  prevent infinite loops which would be
  72.     caused by the repetitive application of certain strategy rules.
  73.     See the section  S T A T E  D E F I N I T I O N  below and see
  74.     section V of the rules.
  75.  
  76.     Unlike the  other  symbolic  transforms, these rules are main-
  77.     tained in a list called InvZTransformRules.   This  was neces-
  78.     sary because Mathematica  reordered these rules in undesirable
  79.     ways if they were coded  as a set of transformation functions.
  80.     Another benefit of  the list form is that it allowed more con-
  81.     trol over how rules are applied.
  82.  
  83.     In general,  this rule base works on expressions with negative
  84.     powers of z.  Great  pains  are  taken  to  convert  any input
  85.     expressions  to  this form accordingly.  This conversion works
  86.     very well for rational z-transform functions.
  87.  
  88.     The  InvZTransform  function  first  calls the interface rules
  89.     (InvZTransformInterfaceRules), which handle default arguments,
  90.     error messages, and multidimensional transforms.   These rules
  91.     simply call  MyInvZTransform  once per dimension of the trans-
  92.     form. This inverse z-transform is biased toward causal systems.
  93.     However,  by specifying the proper region of convergence (ROC),
  94.     anti-causal sequences can be found (see InvZTransform::usage).
  95.  
  96.     The driver for the one-dimensional  rule base is the six argu-
  97.     ment version of MyInvZTransform.  MyInvZTransform[X, z, n, rm,
  98.     rp, s, op] either returns the  inverse transform as a function
  99.     of n or an approximation to  the actual inverse transform.  If
  100.     the rule base can do neither,  which  should  not  happen, the
  101.     interface function cleanup will report any errors.   Arguments
  102.     of MyInvZTransform:
  103.  
  104.     X  z-transform function        rm  Rminus of ROC
  105.     z  z-transform variable        rp  Rplus of ROC
  106.     n  "time"-domain variable    st  local state semaphores
  107.                     op  options
  108.  
  109.     The data tag myinvz is  manipulated by the InvZTransformRules:
  110.     myinvz[X, z, n, rm, rp, st, op, zlist, rmlist, rplist].    The
  111.     last three arguments are only used  by the special multidimen-
  112.     sional rules.
  113.  *)
  114.  
  115. (*  :Functions:    InvZTransform  *)
  116.  
  117.  
  118.  
  119. (*  B E G I N     P A C K A G E  *)
  120.  
  121. BeginPackage[ "SignalProcessing`Digital`InvZTransform`",
  122.           "SignalProcessing`Digital`ZSupport`",
  123.           "SignalProcessing`Support`TransSupport`",
  124.           "SignalProcessing`Support`ROC`",
  125.           "SignalProcessing`Support`SigProc`",
  126.           "SignalProcessing`Support`FilterSupport`",
  127.           "SignalProcessing`Support`SupCode`" ]
  128.  
  129.  
  130. If [ TrueQ[ $VersionNumber >= 2.0 ],
  131.      Off[ General::spell ];
  132.      Off[ General::spell1 ] ];
  133.  
  134.  
  135. (*  U S A G E     I N F O R M A T I O N  *)
  136.  
  137. InvZTransform::usage =
  138.     "InvZTransform[e,z] or InvZTransform[e,z,n] gives the inverse \
  139.     z-transform of e. \
  140.     Here, e can be a function of z, a list, or a z-transform object \
  141.     (see ZTransData). \
  142.     As a list, e should contain three elements:  F[z], rminus, and rplus, \
  143.     such that the region of convergence is rminus < |z| < rplus. \
  144.     If e is a z-transform object, then the inverse z-transform is simply \
  145.     InvZTransform[e] because e contains all necessary information. \
  146.     This routine returns signal processing expressions. \
  147.     Convert them to mathematical formulas by using TheFunction."
  148.  
  149. (*  E N D     U S A G E     I N F O R M A T I O N  *)
  150.  
  151.  
  152. Begin["`Private`"]
  153.  
  154.  
  155. (*  U S E R     I N T E R F A C E  *)
  156.  
  157. (*  InvZ operator  *)
  158. Unprotect[InvZ]
  159. InvZ/: TheFunction[ InvZ[z_, n_][f_] ] := InvZTransform[f, z, n]
  160. Protect[InvZ]
  161.  
  162. (*  InvZTransform  *)
  163. InvZTransform/: Options[InvZTransform] :=
  164.     { Apart -> Rational, Dialogue -> True, Simplify -> True,
  165.       Terms -> 10, TransformLookup -> {} }
  166.  
  167. InvZTransform[x_] :=
  168.     invzdriver[ Options[InvZTransform], x ]
  169. InvZTransform[x_, z_] :=
  170.     invzdriver[ Options[InvZTransform], x, z ]
  171. InvZTransform[x_, z_, n_] :=
  172.     invzdriver[ Options[InvZTransform], x, z, n ]
  173. InvZTransform[x_, z_, n_, op__] :=
  174.     invzdriver[ ToList[op] ~Join~ Options[InvZTransform], x, z, n ]
  175.  
  176. (*  Global variables  *)
  177. completeOptions = {}
  178. dialogueFlag = False
  179. dialogueAllFlag = False
  180. indentationString = ""
  181. simplifyFlag = False
  182.  
  183. nList = {}
  184. zList = {}
  185. multidFlag = False
  186. rmList = {}
  187. rpList = {}
  188.  
  189. (*  validvarQ  *)
  190. validvarQ[ x_Symbol ] := True
  191. validvarQ[ x_List ] := Apply[And, Map[VariableQ, x]]
  192. validvarQ[ x_ ] := False
  193.  
  194. (*  Driver for inverse z-transform interface  *)
  195. invzdriver[op_, x_, z_, n_] :=
  196.     Message[ Transform::badvar, "z", InvZTransform, z ] /;
  197.     ! validvarQ[z]
  198.  
  199. invzdriver[op_, x_, z_, n_] :=
  200.     Message[ Transform::badvar, "discrete-time", InvZTransform, n ] /;
  201.     ! validvarQ[n]
  202.  
  203. invzdriver[op_, args__] :=
  204.     Block [ {},
  205.         dialogue = Replace[Dialogue, op];
  206.  
  207.                 (* Set all global variables right here *)
  208.         completeOptions = op;
  209.         dialogueFlag = TrueQ[dialogue] || SameQ[dialogue, All];
  210.         dialogueAllFlag = SameQ[dialogue, All];
  211.         simplifyFlag = TrueQ[Replace[Simplify, op]];
  212.         nList = zList = {};
  213.         multidFlag = False;
  214.  
  215.                 (* Do the transform  *)
  216.         cleanup [ invztransform[op, args], InformUserQ[op] ] ]
  217.  
  218. invztransform[op_, args__] :=
  219.     Replace [ zinverse[op, args], InvZTransformInterfaceRules ]
  220.  
  221. InvZTransformInterfaceRules = {
  222.     zinverse[op_, e_] :>
  223.         invztransform[ op, e, ZVariables[e] ] /; ZTransformQ[e],
  224.     zinverse[op_, e_] :>
  225.         Message[ InvZTransform::novariables, "z", GetVariables[e] ],
  226.  
  227.     zinverse[op_, e_, z_] :>
  228.         invztransform[op, e, z, DummyVariables[Length[z], Global`n]],
  229.  
  230.     zinverse[op_, e_List, z_, n_] :>
  231.         invztransform[op, MakeZObject[e, z], z, n],
  232.     zinverse[op_, e_, z_Symbol, n_] :>
  233.         MyInvZTransform[e, z, n, op],
  234.     zinverse[op_, e_, z_List, n_] :>
  235.         multiDInvZTransform[e, z, n, op]
  236. }
  237.  
  238. cleanup[ Null, flag_ ] := Null
  239.  
  240. cleanup[ trans_, flag_ ] :=
  241.     Block [    {},
  242.         If [ flag, explain[trans]; Scan[explain, trans, Infinity] ];
  243.         trans ]
  244.  
  245. explain[ myinvz[f_, z_, rest__] ] :=
  246.     Message[ Transform::incomplete, "inverse z-transform", f, z ]
  247.  
  248.  
  249. (*  M E S S A G E S  *)
  250.  
  251. InvZTransform::badterms = "Non-integer terms for series expansion: ``"
  252. InvZTransform::badROC = "Improper region of convergence in ``."
  253.  
  254.  
  255. (*  T R A N S F O R M A T I O N     S T A T E  *)
  256.  
  257.     (*  Listed in order of appearance in the rule base  *)
  258.  
  259. polyfactorfield = 1        (* factor denominator of polynomials         *)
  260. expcheckfield = 2        (* check to see if r^n form             *)
  261. normdenomfield = 3        (* normalize denominator             *)
  262. partialfractionsfield = 4    (* partial fractions on poly w Apart         *)
  263. partialfractionsfield2 = 5    (* partial fractions on poly w MyApart         *)
  264. logformfield = 6        (* factoring inside of a logarithmic form    *)
  265. expformfield = 7        (* factoring inside of an exponential form   *)
  266. expandallfield = 8        (* apply ExpandAll[] to current expression   *)
  267. anticausalfield = 9        (* try anti-causal expansion             *)
  268. cepstrumfield = 10        (* apply rule for complex cepstrum to log    *)
  269. seriesexpansionfield = 11    (* power series expansion of cur. expression *)
  270.  
  271. statevariables = 11        (* number of state variables *)
  272.  
  273. (*
  274.     return the initial transformation state, which is a list of whether or not
  275.     partial fractions, the first denominator normalization, and the second
  276.     denominator initialization has been performed yet on that expression.
  277.  *)
  278.  
  279. initInvZstate[] := Table[True, {statevariables}]
  280. nullInvZstate[] := Table[False, {statevariables}]
  281. anticausalInvZstate[] :=
  282.     Block [    {state},
  283.         state = SetStateField [    nullInvZstate[], logformfield, True ];
  284.         SetStateField [ state, expformfield, True ] ]
  285.  
  286. resetstate[s_] := SetStateField[s, expcheckfield, True]
  287.  
  288.  
  289. (*  S U P P O R T I N G     R O U T I N E S  *)
  290.  
  291. (*    Since the inverse z-transform rule base is biased toward causal    *)
  292. (*  sequences, the leftOrRightSided rules adjust the right-sided    *)
  293. (*  inverse returned by InvZTransform.                    *)
  294. (*    One of the strategies handles the general case.            *)
  295.  
  296. leftsided[p_, n_] :=
  297.     ( p /. { Step[b_. n + c_.] :> - Step[- b n - b + c] } )
  298. rightsided[p_, n_] := p
  299.  
  300. leftOrRightSided[a_, p_, n_, rm_, rp_] := leftsided[p, n]  /; SameQ[Abs[a], rp]
  301. leftOrRightSided[a_, p_, n_, rm_, rp_] := leftsided[p, n]  /; SameQ[Abs[-a], rp]
  302. leftOrRightSided[a_, p_, n_, rm_, rp_] := leftsided[p, n]  /; N[Abs[a] > rp]
  303. leftOrRightSided[a_, p_, n_, rm_, rp_] := rightsided[p, n] /; SameQ[Abs[a], rm]
  304. leftOrRightSided[a_, p_, n_, rm_, rp_] := rightsided[p, n] /; SameQ[Abs[-a], rm]
  305. leftOrRightSided[a_, p_, n_, rm_, rp_] := rightsided[p, n] /; N[Abs[a] < rm]
  306. leftOrRightSided[a_, p_, n_, rm_, rp_] := rightsided[p, n] /; InfinityQ[rp]
  307. leftOrRightSided[a_, p_, n_, rm_, rp_] := rightsided[p, n]
  308.  
  309. (* downsampledRewrite *)
  310. downsampledRewrite[ f_, k_, M_, z_ ] :=
  311.     Block [ {},
  312.         downsampledExpRewrite[f, k, M, z] /.
  313.           ( z^b_. :> takeRealPower[z^b, M] ) ]
  314.  
  315. (* downsampledExpRewrite  *)
  316. downsampledExpRewrite[ f_, k_, M_, z_ ] :=
  317.     Block [ {},
  318.         f /. ( z^a_. Exp[c_ k / M] :> z^a /;
  319.              FreeQ[a, z] && (c == -a M I Pi 2) ) ]
  320.  
  321. (* downsampledQ  *)
  322. downsampledQ[ f_, k_, M_, z_ ] :=
  323.     (! FreeQ[f, z]) && (! FreeQ[f, M]) &&
  324.       (! SameQ[ f, downsampledExpRewrite[f, k, M, z] ])
  325.  
  326. (*  fixUp -- strip redundant options and remove TransformLookup option  *)
  327. fixUp[ op_ ] :=
  328.     { Dialogue -> Replace[Dialogue, op],
  329.       Simplify -> Replace[Simplify, op],
  330.       Terms -> Replace[Terms, op] }
  331.  
  332. (*  mDresamplingCheck  *)
  333. mDresamplingCheck[Upsample, x_, z_] :=
  334.     Block [ {cond, upmatrixtrans, zvars},
  335.         cond = ! MyFreeQ[ GetAllFactors[x, z], zList ];
  336.         If [ cond,
  337.              Needs[ "SignalProcessing`Digital`MDInvZTransform`" ];
  338.              cond = upsampleSystemQ[x, zList] ];
  339.         cond ]
  340.  
  341. (*  multByExponential  *)
  342. multByExponentialQ[f_, z_] :=
  343.     Block [    {scale},
  344.         scale = ScalingFactor[f, z];
  345.         ( ! SameQ[scale, 1] ) && ( ! SameQ[scale, 0] ) ]
  346.  
  347. (*  multiDInvZTransform  *)
  348. multiDInvZTransform[e_, z_, n_, op_] :=
  349.     Block [    {invztrans, oldflag, oldnList, oldrmlist, oldrplist, oldzList},
  350.         {oldzList, oldnList, oldFlag, oldrmlist, oldrplist} =
  351.             {zList, nList, multidFlag, rmList, rpList};
  352.         zList = z;
  353.         nList = n;
  354.         multidFlag = True;
  355.         If [ ZTransformQ[e],
  356.              {rmList, rpList} = { GetRMinus[e], GetRPlus[e] },
  357.              {rmList, rpList} = { {}, {} } ];
  358.         invztrans = MultiDInvTransform [ e, z, n, op, ZTransformQ,
  359.                          MyInvZTransform, MakeZObject,
  360.                          Global`n ];
  361.         {zList, nList, multidFlag, rmList, rpList} =
  362.             {oldzList, oldnList, oldFlag, oldrmlist, oldrplist};
  363.         invztrans ]
  364.  
  365. (*  normalizeDenominator  *)
  366. normalizeDenominator[ratpoly_, z_] :=
  367.     Block [    {adjdenom, adjnumer, const, denom, maxexp, numer, term},
  368.         numer = Numerator[ratpoly];
  369.         denom = Expand[Denominator[ratpoly]];
  370.         maxexp = Exponent[denom, z];
  371.         term = z^(-maxexp);
  372.         adjdenom = Distribute[term denom];
  373.         const = ConstantTerm[adjdenom, z];
  374.         adjdenom = Distribute[ (1/const) adjdenom ];
  375.         adjnumer = If [ MixedPolynomialQ[numer, z],
  376.                 Distribute[ term Expand[numer] ],
  377.                 term numer ];
  378.         1/const adjnumer / adjdenom ]
  379.  
  380. (*  partialFractionsDialogue  *)
  381. partialFractionsDialogue[p_, partfrac_, options_, apart_] :=
  382.     Block [ {},
  383.         If [ InformUserQ[options],
  384.              Print["( after the term"];
  385.              Print["  ", p];
  386.              Print["  has been broken up into its partial " ];
  387.              Print["  fractions representation using ", apart];
  388.              Print["  ", partfrac, " )" ] ] ]
  389.  
  390. (*  partialFractionsQ  *)
  391. partialFractionsQ[p_, z_] :=
  392.     MixedPolynomialQ[ Expand[Denominator[p]], z ] &&
  393.       ! FreeQ[ Denominator[p], z ]
  394.  
  395. (*  rootROC  *)
  396. SetAttributes[rootROC, {Listable}]
  397. rootROC[x_, y_] := takeRealPower[x, y]
  398.  
  399. (*  seriesDialogue  *)
  400. seriesDialogue[x_, expansion_, terms_, options_] :=
  401.     Block [    {},
  402.         If [ InformUserQ[options],
  403.              Print["( after breaking up the expression"];
  404.              Print["  ", x];
  405.              Print["  into its series representation"];
  406.              Print["  ", expansion, " )"] ] ]
  407.  
  408. (*  seriesExpansion  *)
  409. seriesExpansion[ True, x_, z_, power_ ] := Series[x, {z, 0, power}]
  410.  
  411. seriesExpansion[ False, x_, z_, power_ ] :=
  412.     Block [    {negx, expansion},
  413.         negx = x /. z -> z^-1;
  414.         expansion = Series[negx, {z, 0, power}];
  415.         expansion /. z -> z^-1 ]
  416.  
  417. (*  seriesStrategy  *)
  418. seriesStrategy[ x_, z_, n_, rm_, rp_, s_, op_, rest___ ] :=
  419.     Block [ {expansion, exponents, posexpandflag, state, terms},
  420.         terms = Replace[Terms, op];
  421.         state = SetStateField[s, seriesexpansionfield, False];
  422.  
  423.         (*  If the Terms option is disabled, then take no action.  *)
  424.         If [ ! terms,
  425.              Return[ myinvz[x, z, n, rm, rp, state, op, rest] ] ];
  426.  
  427.         If [ ! IntegerQ[terms],
  428.              Message[ InvZTransform::badterms, terms ];
  429.              myinvz[x, z, n, rm, rp, state, op, rest] ];
  430.  
  431.         (*  See if we should expand in pos. or neg. powers    *)
  432.         exponents = GetAllExponents[x, z];
  433.         posexpandflag = If [ Apply[And, Thread[exponents >= 0]],
  434.                      True, False, False ];
  435.  
  436.         (*  First expansion  *)
  437.         expansion = seriesExpansion[posexpandflag, x, z, terms - 1];
  438.  
  439.         (*  Second expansion if first failed  *)
  440.         If [ ! SameQ[Head[expansion], SeriesData],
  441.              expansion =
  442.                seriesExpansion[Not[posexpandflag], x, z, terms - 1] ];
  443.  
  444.         (*  If this fails, call MyInvZTransform again        *)
  445.         (*  but with series expansion disabled            *)
  446.         If [ ! SameQ[Head[expansion], SeriesData],
  447.              myinvz[ x, z, n, rm, rp, state, op, rest ],
  448.              seriesDialogue[x, expansion, terms, op];
  449.              state = nullInvZstate[];
  450.              Map [ myinvz[#1, z, n, rm, rp, state, op, rest]&, 
  451.                Normal[expansion] ] ] ]
  452.  
  453. (*  takeRealPower --  assume that both arguments are real-valued  *)
  454. SetAttributes[takeRealPower, {Listable}]
  455. takeRealPower[ 0, r_ ] := 0
  456. takeRealPower[ Infinity, r_ ] := Infinity
  457. takeRealPower[ x_^b_, r_ ] := x^(b r) /; IntegerQ[b r]
  458. takeRealPower[ x_, r_ ] := x^r
  459.  
  460. (*  upsampledZtransform  *)
  461. upsampledZtransform[ f_, z_, n_, rm_, rp_, s_, op_, rest___] :=
  462.     Block [    {newf, newrm, newrp, upindex},
  463.         upindex = UpsampleFactor[f, z];
  464.         If [ IntegerQ[upindex],
  465.              upindex = Abs[upindex],
  466.              If [ SameQ[Replace[Dialogue, op], All],
  467.               upindex = Replace[upindex, -k_ :> k];
  468.                   Print[ "assuming ", upindex, " is an integer" ] ] ];
  469.         newf = f /. ( z^a_. :> z^(a/upindex) );
  470.         newrm = rootROC[rm, upindex];     (* adjust ROC *)
  471.         newrp = rootROC[rp, upindex];
  472.         Upsample[upindex, n]
  473.             [myinvz[newf, z, n, newrm, newrp, s, op, rest]] ]
  474.  
  475.  
  476. (*  B E G I N     R U L E     B A S E  *)
  477.  
  478. ROCrewriterules = {
  479.     b_?Positive / Abs[a_] :> Abs[b / a],
  480.     b_?Negative / Abs[a_] :> - Abs[- b / a]
  481. }
  482.  
  483. (*  MyInvZByPass--- called by inverse DTFT and interfaces to MyInvZTransform  *)
  484. MyInvZByPass[ e_, z_, n_, zlist_, nlist_, op_ ] :=
  485.     MyInvZTransform[ e, z, n, op ]
  486.  
  487. (*  MyInvZTransform  *)
  488. MyInvZTransform[ e_, z_, n_, op_ ] :=
  489.     Block [ {rminus, rplus, valid},
  490.         rminus = SPSimplify[ GetRMinus[e] ];
  491.         rplus = SPSimplify[ GetRPlus[e] ] /. ROCrewriterules;
  492.         valid = Apply[And, Thread[ToList[rminus] < ToList[rplus]]];
  493.         If [ TrueQ[! valid],
  494.              Message[ InvZTransform::badROC, e ],
  495.              MyInvZTransform[ TheFunction[e], z, n, rminus,
  496.                            rplus, initInvZstate[], op ] ] ] /;
  497.     ZTransformQ[e]
  498.  
  499. MyInvZTransform[ e_, z_, n_, op_ ] :=
  500.     MyInvZTransform[ e, z, n, 0, Infinity, initInvZstate[], op ]
  501.  
  502.  
  503. (*   Driver for one-dimensional rule base            *)
  504. (*   First, convert all z /(z - a) forms to 1/(1 - a z^-1)    *) 
  505. (*    since the rule base favors terms with z^-1 terms    *)
  506. (*   Loop until the expression to be inverted does not change.    *)
  507. MyInvZTransform[ e_, z_, n_, rm_, rp_, st_, options_ ] :=
  508.     Block [ {laste, newe, newinvzrules, nvars, op, trace},
  509.  
  510.                 (* generate the z-transform rules *)
  511.         newinvzrules = TransformFixUp[ z, n, options, myinvz, False,
  512.                            InvZTransform, Null, Null ];
  513.         InvZTransformRules = Join[ newinvzrules,
  514.                        OriginalInvZTransformRules ];
  515.  
  516.         (* strip redundant options and set dialogue level *)
  517.         op = fixUp[options];
  518.         trace = SameQ[ Replace[Dialogue, op], All ];
  519.  
  520.         (* take the 1-D inverse transform *)
  521.         newe = myinvz[ MapAll[convert[#, z]&, e],
  522.                    z, n, rm, rp, st, op, zList, rmList, rpList ];
  523.         While [ ! SameQ[ laste, newe ],
  524.             If [ trace, Print[ newe ]; Print[ "which becomes" ] ];
  525.             laste = newe;
  526.             newe = laste /. ( myinvz[a__] :>
  527.                  Replace[myinvz[a], InvZTransformRules] ) ];
  528.  
  529.         (* fix up the result *)
  530.         newe = laste /. postinvzrules;
  531.         While [ ! SameQ[ laste, newe ],
  532.             If [ trace, Print[ newe ]; Print[ "which becomes" ] ];
  533.             laste = newe;
  534.             newe = laste /. postinvzrules ];
  535.         If [ trace, Print[ newe ] ];
  536.  
  537.         (* simplify the resulting expression if requested *)
  538.         If [ Replace[Simplify, op],
  539.              nvars = If [ EmptyQ[nList], {n}, nList ];
  540.              newe = SPSimplify[ newe,
  541.                     Dialogue -> trace,
  542.                     Variables -> nvars ] ];
  543.  
  544.         newe ]
  545.  
  546. (*  normalize as much of the expression as possible first  *)
  547. convert[ z_^k_ ( z_ + a_ )^j_, z_ ] := ( 1 + a z^-1 )^j  /; j == -k
  548. convert[ z_ / ( z_ + a_ ), z_ ] := 1 / ( 1 + a z^-1 )
  549. convert[ x_, z_ ] := x
  550.  
  551. postinvzrules = {
  552.     shift[a_. t_, n_, n0_] :>
  553.         If [ Count[{t}, f_[p___][b___], Infinity] > 0,
  554.              a Shift[n0, n][t],            (* operators present *)
  555.              a (t /. n -> n - n0) ] /;        (* no operators *)
  556.         FreeQ[a, n],
  557.  
  558.     (*  These two rules are needed to support a series    *)
  559.     (*  expansion which returns a constant term plus    *)
  560.     (*  terms of the form  approx z^r  .            *)
  561.     myinvz[ b_, z_, n_, rest___ ] :>
  562.         b Impulse[n] /;
  563.         FreeQ[b, z],           (* Converges all ROC *)
  564.  
  565.     myinvz[ a_. z_^r_., z_, n_, rest___ ] :>
  566.         a Impulse[n + r] /;
  567.         FreeQ[{a,r}, z],       (* Converges all ROC *)
  568.  
  569.     (*  Attempt a series expansion about z = 0        *)
  570.     (*  This is the strategy when all else has failed.    *)
  571.     (*  The user can disable this rule (Terms -> False).    *)
  572.     (*  Please keep this as the last rule.            *)
  573.     myinvz[x_, z_, n_, rm_, rp_, s_, op_, rest___] :>
  574.         seriesStrategy[x, z, n, rm, rp, s, op, rest ] /;
  575.         GetStateField[s, seriesexpansionfield]
  576. }
  577.  
  578. Format[ myinvz[ x_, z_, n_, rest___ ] ] := 
  579.     SequenceForm[ ColumnForm[{"Z" Superscript[-1],
  580.                   "  " ~StringJoin~ ToString[z]}],
  581.               { x } ]
  582. Format[ shift[t_, n_, n0_] ] := SequenceForm[ {t}, Subscript[n -> n + n0] ]
  583.  
  584.  
  585. (*  Most of the rest of the file is the list of inverse z-transform rules  *)
  586.  
  587. InvZTransformRules = {}
  588.  
  589. OriginalInvZTransformRules = {
  590.  
  591.  
  592.   (*  M U L T I D I M E N S I O N A L     H O O K S  *)
  593.  
  594.   myinvz[ f_, z_, n_, rm_, rp_, s_, op_, rest___ ] :>
  595.     upsampleSetUp[ f, z, n, rm, rp, s, op, rest, nList ] /;
  596.     multidFlag && mDresamplingCheck[Upsample, f, z],
  597.  
  598.   myinvz[ newROC[ f_, args___ ], z_, n_, rest___ ] :>
  599.     newROC[ myinvz[ f, z, n, rest ], args ],
  600.  
  601.  
  602.  
  603.  
  604.   (*  I.  R A T I O N A L     T R A N S F O R M S        *)
  605.  
  606.  
  607.   (*      A.  Simple transform pairs                *)
  608.   (*          1.  Constant (handles b=0 case as well)        *)
  609.   myinvz[ b_, z_, n_, rest___ ] :>
  610.     b Impulse[n] /;
  611.     FreeQ[b, z],                    (* Converges all ROC *)
  612.  
  613.   (*              Redundant rule --  placed here for speed    *)
  614.   myinvz[ a_. z_^r_., z_, n_, rest___ ] :>
  615.     a Impulse[n + r] /;
  616.     FreeQ[{a,r}, z],                (* Converges all ROC *)
  617.  
  618.   (*          2.  Discrete-time sine functions            *)
  619.   myinvz[ 1 / (1 - 2 Cos[w0_] z_^-1 + z_^-2), z_, n_, rm_, rp_, rest___ ] :>
  620.     leftOrRightSided[1, Sin[w0 + w0 n] Step[n] / Sin[w0], n, rm, rp] /;
  621.     FreeQ[w0, z],
  622.  
  623.   myinvz[ 1 / (1 + c_. z_^-1 + z_^-2), z_, n_, rm_, rp_, rest___ ] :>
  624.     Block [    {sinw0, w0},
  625.         w0 = ArcCos[- c / 2];
  626.         sinw0 = Sqrt[1 - c^2/4];
  627.         leftOrRightSided[ 1, Sin[w0 n + w0] Step[n] / sinw0,
  628.                   n, rm, rp ] ] /;
  629.     N[-2 < c < 2],
  630.  
  631.   (*          3.  Discrete-time cosine functions        *)
  632.  
  633.   myinvz[ (1 - Cos[w0_] z_^-1) / (1 - 2 Cos[w0_] z_^-1 + z_^-2), z_, n_,
  634.       rm_, rp_, rest___ ] :>
  635.     leftOrRightSided[1, Cos[w0 n] Step[n], n, rm, rp ] /;
  636.     FreeQ[w0, z],
  637.  
  638.   (*          4.  Discrete-time hyperbolic functions        *)
  639.  
  640.   myinvz[ (Sinh[b_] + Sinh[c_] z_^-1) / (1 - 2 Cosh[w_] z_^-1 + z_^-2), z_, n_,
  641.       rm_, rp_, rest___ ] :>
  642.     leftOrRightSided[1, Sinh[b + w n] Step[n], n, rm, rp] /;
  643.     FreeQ[{w,b,c},z] && ( c == w - b ),
  644.  
  645.   myinvz[ (Cosh[b_] + Cosh[c_] z_^-1) / (1 - 2 Cosh[w_] z_^-1 + z_^-2), z_, n_,
  646.       rm_, rp_, rest___ ] :>
  647.     leftOrRightSided[1, Cosh[b + w n] Step[n], n, rm, rp] /;
  648.     FreeQ[{w,b,c},z] && ((c == w - b) || (c == b - w)),
  649.  
  650.  
  651.   (*      B.  First order sections                    *)
  652.   (*          most covered by multiplication-by-exponential property    *)
  653.   (*          1/(1+az^-1) is more general than 1/(1-az^-1)        *)
  654.  
  655.   (*          1.  Discrete-time pulse has no poles            *)
  656.   myinvz[ ( 1 + b_. z_^L_ ) / ( 1 + a_. z_^-1 ), z_, n_, rm_, rp_, rest___ ] :>
  657.     Block [ {},
  658.         Assuming[ L < 0, dialogueAllFlag ];
  659.         If [ dialogueAllFlag, Print[ "where ", L, " is an integer" ] ];
  660.         (-a)^n Pulse[-L, n] ] /;
  661.     FreeQ[{a,b,L}, z] && Implies[ NumberQ[L], IntegerQ[L] && L < 0 ] &&
  662.       ( Equal[a, b, 1] || SameQ[-b, (-a)^(-L)] ),
  663.  
  664.   (*          2.  General exponential times a step (right-sided pair)    *)
  665.   myinvz[ 1 / ( 1 + a_. z_^-1 ), z_, n_, rm_, rp_, rest___ ] :>
  666.     leftOrRightSided[a, (-a)^n Step[n], n, rm, rp] /;
  667.     FreeQ[a,z],
  668.  
  669.   (*          3.  General first-order section                *)
  670.   myinvz[ a_. / (c_ + d_. z_^-1), z_, n_, rm_, rp_, rest___ ] :>
  671.     leftOrRightSided[d/c, a/c (-d/c)^n Step[n], n, rm, rp] /;
  672.     FreeQ[{a,c,d}, z],
  673.  
  674.   myinvz[ (a_. + b_. z_^-1) / (c_ + d_. z_^-1), z_, n_, rm_, rp_, rest___ ] :>
  675.     leftOrRightSided[d/c,
  676.              a/c (-d/c)^n Step[n] + b/c (-d/c)^(n-1) Step[n-1],
  677.              n, rm, rp] /;
  678.     FreeQ[{a,b,c,d}, z],
  679.  
  680.  
  681.   (*      C.  Second order sections                    *)
  682.  
  683.   myinvz[ z_^-1 / ( 1 + b_. z_^-1 + a_. z_^-2 ), z_, n_, rm_, rp_, rest___ ] :>
  684.     Block [    {costheta, fun, rho, theta},
  685.         rho = Sqrt[a];
  686.         costheta = -b / ( 2 rho );
  687.         theta = ArcCos[costheta];
  688.         fun = ( 1 / ( rho Sin[theta] ) ) rho^n Sin[n theta] Step[n];
  689.         leftOrRightSided[Sqrt[a], fun, n, rm, rp] ] /;
  690.     N[(a > 0) && (4 a > b^2)],
  691.  
  692.   myinvz[ (1 + b_. z_^-1) / (1 + c_. z_^-1 + a_. z_^-2), z_, n_,
  693.       rm_, rp_, rest___ ] :>
  694.     Block [    {costheta, fun, rho, theta},
  695.         rho = Sqrt[a];
  696.         costheta = - b / rho;
  697.         theta = ArcCos[costheta];
  698.         fun = rho^n Cos[n theta] Step[n];
  699.         leftOrRightSided[Sqrt[a], fun, n, rm, rp] ] /;
  700.     N[ (c == 2 b) && (0 < a < b^2) ],
  701.  
  702.   myinvz[ 1 / ( 1 + b_. z_^-1 + a_. z_^-2 ), z_, n_, rm_, rp_, rest___ ] :>
  703.     Block [ {costheta, rho, sintheta, theta},
  704.         rho = Sqrt[a];
  705.         costheta = - b / ( 2 rho );
  706.         theta = ArcCos[costheta];
  707.         sintheta = Sin[theta];
  708.         fun = rho^n ( Cos[n theta] + Cot[theta] Sin[n theta] ) Step[n];
  709.         leftOrRightSided[Sqrt[a], fun, n, rm, rp] ] /;
  710.     N[ (a > 0) && (4 a > b^2) ],
  711.  
  712.   myinvz[ 1 / ( 1 + a_. z^-2 ), z_, n_, rm_, rp_, s_, op_, rest___ ] :>
  713.     Block [ {},
  714.         Assuming[ Negative[a], dialogueAllFlag ];
  715.         leftOrRightSided[ Sqrt[-a],
  716.                   ( Sqrt[-a]^n + (- Sqrt[-a])^n ) Step[n] / 2,
  717.                   n, rm, rp ] ] /;
  718.     FreeQ[a,z] && Implies[ NumberQ[N[a]], N[a < 0] ],
  719.  
  720.   myinvz[ 1 / ( 1 + a_. z_^-2 ), z_, n_, rm_, rp_, s_, op_, rest___ ] :>
  721.     Block [ {},
  722.         Assuming[ Positive[a], dialogueAllFlag ];
  723.         leftOrRightSided[ Sqrt[a],
  724.                   Sqrt[a]^n Cos[n Pi / 2] Step[n],
  725.                   n, rm, rp ] ] /;
  726.     FreeQ[a,z] && Implies[ NumberQ[N[a]], N[a > 0] ],
  727.  
  728.   myinvz[ 1 / ( 1 + a_. z_^-1 )^2, z_, n_, rm_, rp_, rest___ ] :>
  729.     leftOrRightSided[a, (n + 1) (-a)^n Step[n], n, rm, rp] /;
  730.     FreeQ[a,z],
  731.  
  732.  
  733.   (*      D.  Higher order sections                    *)
  734.   (*          Note that k can be real-valued.                *)
  735.  
  736.   (*          1.  General denominator section to an integer power    *)
  737.   myinvz[ c_. (1 + b_. z_^-1)^k_, z_, n_, rm_, rp_, rest___ ] :>
  738.     leftOrRightSided[b, c Binomial[n-k-1,-k-1] (-b)^n Step[n], n, rm, rp] /;
  739.     FreeQ[{b,c,k}, z] && ( k < -1 ),
  740.  
  741.   (*          2.  General denominator section to a symbolic power    *)
  742.   (*          Step is different because we don't know what k is.    *)
  743.   myinvz[ c_. / (1 + b_. z_^-1)^k_., z_, n_, rm_, rp_, rest___ ] :>
  744.     leftOrRightSided[b, c Binomial[n+k-1,k-1] (-b)^n Step[n+k-1],
  745.              n, rm, rp] /;
  746.     FreeQ[{b,c,k}, z],
  747.  
  748.   (*          3.  General numerator section to a symbolic power        *)
  749.   myinvz[ c_. (a_ + b_. z_^-1)^k_, z_, n_, rm_, rp_, s_, op_, rest___ ] :>
  750.     Block [ {},
  751.         Assuming[ Positive[k], dialogueAllFlag ];
  752.         c leftOrRightSided[ -b/a, Binomial[k,n] b^n a^(k-n) Step[n],
  753.                     n, rm, rp ] ] /;
  754.     FreeQ[{a,b,c,k}, z] && Implies[NumberQ[k], N[k > 0]],
  755.  
  756.   (*          4.  Reduce expressions in z^k to expressions in z.    *)
  757.   (*          This applies to expressions other than polys in z.    *)
  758.   myinvz[ f_, z_, n_, rest___ ] :>
  759.     upsampledZtransform[ f, z, n, rest ] /;
  760.     ZUpsampledQ[f, z],
  761.  
  762.  
  763.   (*      E.  Factor rational polynomials and make the denominator    *)
  764.   (*        only contain negative powers of z.            *)
  765.  
  766.   (*          1.  Factor rational polynomials when denominators have    *)
  767.   (*          positive powers of z (ensured by PolynomialQ test).    *)
  768.   myinvz[ p_, z_, n_, rm_, rp_, s_, rest___ ] :>
  769.     myinvz [ ( Expand[Numerator[p]] ) / ( Factor[Denominator[p]] ),
  770.          z, n, rm, rp,
  771.          SetStateField[s, polyfactorfield, False], rest ] /;
  772.     GetStateField[s, polyfactorfield] && PolynomialQ[Denominator[p], z] &&
  773.       RationalFunctionQ[p, z],
  774.  
  775.   (*          2.  Make factored denominators have neg. powers of z    *)
  776.   (*          a.  For numeric k, numeric j                *)
  777.   myinvz[ ((a_ + b_. z_^j_.)^k_) x_., z_, n_, rest___ ] :>
  778.     b^k myinvz[ Distribute[Numerator[x] z^(j k)] ( 1 + a z^(-j) / b )^k /
  779.             Denominator[x], z, n, rest ] /;
  780.     FreeQ[{a,b}, z] && IntegerQ[k] && (k < 0) && IntegerQ[j] && (j > 0),
  781.  
  782.   (*          b.  For symbolic k, numeric j                *)
  783.   myinvz[ x_. / ((a_ + b_. z_^j_.)^k_), z_, n_, rest___ ] :>
  784.     b^-k myinvz[ Distribute[Numerator[x] z^-(j k)] / (1 + a z^(-j) / b)^k /
  785.              Denominator[x], z, n, rest ] /;
  786.     FreeQ[{a,b,k}, z] && IntegerQ[j] && (j > 0),
  787.  
  788.   (*          3.  Normalize factored denominators.            *)
  789.   (*          a.  For numeric k, numeric j                *)
  790.   myinvz[ ((a_ + b_. z_^j_)^k_) x_., z_, n_, rest___ ] :>
  791.     a^k myinvz[ x ( 1 + b z^j / a ) ^ k, z, n, rest ] /;
  792.     FreeQ[{a,b}, z] && IntegerQ[k] && (k < 0) && IntegerQ[j] &&
  793.       (j < 0) && ! SameQ[a, 1],
  794.  
  795.   (*          b.  For symbolic k, numeric j                *)
  796.   myinvz[ x_. / ((a_ + b_. z_^j_)^k_), z_, n_, rest___ ] :>
  797.     a^-k myinvz[ x / ( 1 + b z^j / a ) ^ k, z, n, rest ] /;
  798.     FreeQ[{a,b,k}, z] && IntegerQ[j] && (j < 0) && ! SameQ[a, 1],
  799.  
  800.  
  801.  
  802.  
  803.   (*  II. N O N - R A T I O N A L     T R A N S F O R M S  *)
  804.  
  805.  
  806.   (*      A.  Handle Exp[affine] forms; covers Sinh[affine]    *)
  807.   (*          Cosh[affine] forms converges for any ROC        *) 
  808.   myinvz[ c_^(b_. + a_. z_^-1), z_, n_, rm_, rp_, rest___ ] :>
  809.     leftOrRightSided[ a Log[c], c^b a^n (Log[c])^n Step[n] / n!,
  810.               n, rm, rp ]  /;
  811.     FreeQ[{a,b,c}, z],
  812.  
  813.  
  814.   (*      B.  Handle Log[affine] forms [O&S, 55]        *)
  815.   myinvz[ Log[a_ + b_. z_^-1], z_, n_, rm_, rp_, rest___ ] :>
  816.     Log[a] Impulse[n] +
  817.       leftOrRightSided[-b/a, -(-b/a)^n Step[n-1] / n, n, rm, rp] /;
  818.     FreeQ[{a,b}, z],
  819.  
  820.  
  821.   (*      C.  Family of causal IIR filters            *)
  822.   myinvz[ ( 1 + z_^-1 )^r_, z_, n_, rm_, rp_, s_, op_, rest___ ] :>
  823.     Block [ {},
  824.         Assuming[ Positive[r + 1], dialogueAllFlag ];
  825.         Gamma[1 + r] Step[n] / ( Gamma[1 + n] Gamma[1 + r - n] ) ] /;
  826.     FreeQ[r,z] && Implies[ NumberQ[N[r]], N[r > -1] ],
  827.  
  828.  
  829.   (*      D.  Other unusual forms.                *)
  830.   myinvz[ ArcTan[z_^-1], z_, n_, rest___ ] :>
  831.     Sin[n Pi / 2] Step[n - 1] / n,
  832.  
  833.   myinvz[ Sqrt[z_] ArcTan[Sqrt[z_^-1]], z_, n_, rest___ ] :>
  834.     Step[n] / ( 2 n + 1 ),
  835.  
  836.   myinvz[ Sqrt[z_] ArcTan[z_^-1], z_, n_, rest___ ] :>
  837.     Step[n] / ( 2 n + 1 ),
  838.  
  839.   myinvz[ Cosh[ Sqrt[z_^-1] ], z_, n_, rest___ ] :>
  840.     Step[n] / ( 2 n ) !,
  841.  
  842.   myinvz[ Sinh[ Sqrt[z_^-1] ], z_, n_, rest___ ] :>
  843.     Step[n] / ( 2 n + 1 ) !,
  844.  
  845.   myinvz[ Sqrt[ 1 / ( 1 + a_. z_^-1 ) ], z_, n_, rest___ ] :>
  846.     (-a)^n (2 n)! Step[n] / ((2^n n!) ^ 2)  /;
  847.     FreeQ[a, z],
  848.  
  849.  
  850.  
  851.  
  852.   (* III. Z - T R A N S F O R M     P R O P E R T I E S            *)
  853.  
  854.   (*  Omitted these properties:    1.  multiplication by n^m        *)
  855.   (*                2.  multiplication by Cos[bn]        *)
  856.   (*                3.  multiplication by Sin[bn]        *)
  857.   (*                4.  multiplication by x y        *)
  858.   (*                                    *)
  859.   (*  This section is similar to the z-transform Properties section    *)
  860.   (*  of the forward z-transform rules.                    *)
  861.  
  862.  
  863.   (*      A.  Downsampled sequence.                    *)
  864.   (*          Had to put it before mult. by exponential sequence    *)
  865.   myinvz[ Summation[k_, 0, Mminus1_, 1][f_], z_, n_, rm_, rp_, rest___ ] :>
  866.     Block [    {m, newf, newrm, newrp},
  867.         m = Mminus1 + 1;
  868.         newf = downsampledRewrite[f, k, m, z];
  869.         newrm = rootROC[rm, 1/m];
  870.         newrp = rootROC[rp, 1/m];
  871.         m Downsample[m, n][myinvz[newf, z, n, newrm, newrp, rest]] ] /;
  872.     downsampledQ[f, k, Mminus1 + 1, z],
  873.  
  874.  
  875.   (*      B.  Multiplication by exponential sequence (exp)^n        *)
  876.   (*          Check this possibility out before partial fractions    *)
  877.   (*          Moved from rational section because it interfered        *)
  878.   (*          with the Log lookup.  This rule should not interfere    *)
  879.   (*          with the downsampling rule which requires a symbolic    *)
  880.   (*          summation to be present.                    *)
  881.   myinvz[ f_, z_, n_, rm_, rp_, s_, rest___ ] :>
  882.     Block [    {exp, fun, newf, newrm, newrp, state},
  883.         exp = 1 / ScalingFactor[f, z];
  884.         newf = f /. z -> z exp;
  885.         newrm = Abs[rm / exp];
  886.         newrp = If [ InfinityQ[rp], Infinity, Abs[rp / exp] ];
  887.         state = SetStateField[s, expcheckfield, False];
  888.         (exp)^n myinvz[newf, z, n, newrm, newrp, state, rest] ] /;
  889.     GetStateField[s, expcheckfield] && FreeQ[f, Summation] &&
  890.     multByExponentialQ[f, z],
  891.  
  892.  
  893.   (*      C.  Homogeneity --  pick off constants            *)
  894.   myinvz[ c_ x_, z_, n_, rest___ ] :>
  895.     c myinvz[ x, z, n, rest ] /;
  896.     FreeQ[c, z],
  897.  
  898.  
  899.   (*      D.  Delay or shift                        *)
  900.   myinvz[ x_ z_^m_., z_, n_, rest___ ] :>
  901.     shift[ myinvz[x, z, n, rest], n, -m ] /;
  902.     FreeQ[m, z] && Implies[ NumberQ[m], RationalQ[m] ],
  903.  
  904.  
  905.   (*      E.  Simple Additivity                        *)
  906.   myinvz[ x_ + y_, z_, n_, rm_, rp_, s_, rest___ ] :>
  907.     myinvz[x, z, n, rm, rp, resetstate[s], rest] +
  908.     myinvz[y, z, n, rm, rp, resetstate[s], rest],
  909.  
  910.  
  911.   (*      F.  Partial fractions decomposition                *)
  912.   (*          partial fractions decomposition for denominators that    *)
  913.   (*          are polynomials in of z; once partial fractions de-    *)
  914.   (*          composition has been applied to an expression, it is    *)
  915.   (*          never applied to any of the resulting sub-expressions.    *)
  916.   (*          This is controlled by the partial fractions semaphore/    *)
  917.   (*          flag in the local state; at this point, the denominator    *)
  918.   (*          should be a function of z^-1, but we will do partial    *)
  919.   (*          fractions on a denominator that is a "two-sided" poly-    *)
  920.   (*          nomial (one in both positive and negative    powers of z).    *)
  921.  
  922.   (*          1.  Normalize the denominator                *)
  923.   myinvz[ p_, z_, n_, rm_, rp_, s_, rest___ ] :>
  924.     myinvz[ normalizeDenominator[p, z], z, n, rm, rp,
  925.         SetStateField[s, normdenomfield, False], rest ] /;
  926.     GetStateField[s, normdenomfield] && partialFractionsQ[p, z],
  927.  
  928.   (*          2.  Use Mathematica's Apart primitive            *)
  929.   myinvz[ p_, z_, n_, rm_, rp_, s_, rest___ ] :>
  930.     Block [ {failure, partfrac, pnorm, pfnorm, state},
  931.         partfrac = KeepNormalized[p, Apart, z, z];
  932.         state = SetStateField[s, partialfractionsfield, False];
  933.         failure = partialFractionsQ[partfrac, z] &&
  934.                   ( SameQ[partfrac, p] ||
  935.                 SameQ[normalizeDenominator[p, z],
  936.                   normalizeDenominator[partfrac, z]] );
  937.         If [ failure,
  938.              partfrac = p,
  939.              partialFractionsDialogue[p, partfrac, op, Apart] ];
  940.         myinvz [ partfrac, z, n, rm, rp, state, rest ] ] /;
  941.     GetStateField[s, partialfractionsfield] && partialFractionsQ[p, z],
  942.  
  943.   (*           3.  Do the partial fractions that Apart cannot do    *)
  944.   myinvz[ p_, z_, n_, rm_, rp_, s_, op_, rest___ ] :>
  945.     Block [ {filter, partfrac, state},
  946.         filter = If [ SameQ[Replace[Apart, op], All], N, Identity ];
  947.         partfrac = KeepNormalized[p, MyApart, z, z, filter];
  948.         state = SetStateField[s, partialfractionsfield, False];
  949.         state = SetStateField[state, partialfractionsfield2, False];
  950.         If [ ! SameQ[partfrac, p],
  951.              partialFractionsDialogue[p, partfrac, op, MyApart] ];
  952.         myinvz [ partfrac, z, n, rm, rp, state, op, rest ] ] /;
  953.     GetStateField[s, partialfractionsfield2] && partialFractionsQ[p, z],
  954.  
  955.  
  956.   (*      G.  Additivity of Numerator                    *)
  957.   myinvz[ (x_+y_)/c_, z_, n_, rm_, rp_, s_, rest___ ] :>
  958.     myinvz[x/c, z, n, rm, rp, resetstate[s], rest] +
  959.     myinvz[y/c, z, n, rm, rp, resetstate[s], rest],
  960.  
  961.  
  962.   (*      H.  Logarithm properties                     *)
  963.   myinvz[ Log[a_. / b_], z_, n_, rest___ ] :>
  964.     myinvz[Log[a] - Log[b], z, n, rest],
  965.  
  966.   myinvz[ Log[a_ b_], z_, n_, rest___ ] :>
  967.     myinvz[Log[a] + Log[b], z, n, rest],
  968.  
  969.  
  970.   (*      I.  Derivative                        *)
  971.   myinvz[ Derivative[m_][f_][z_], z_, n_, rest___ ] :>
  972.     shift[ ZPolynomial[m, n] myinvz[f[z], z, n, rest], n, -m ] /;
  973.     FreeQ[m, z],
  974.  
  975.  
  976.  
  977.  
  978.   (* IV.  S I G N A L     P R O C E S S I N G     S T R U C T U R E S    *)
  979.  
  980.  
  981.   (*        -.  An operator independent of n---  take z-transform of f    *)
  982.   myinvz[ operator_[params__][f_], z_, n_, rest___ ] :>
  983.     operator[params][ myinvz[f, z, n, rest] ] /;
  984.     FreeQ[{params}, z],
  985.  
  986.  
  987.   (*      A.  Difference equations                    *)
  988.   myinvz[ x_ (1 - z_^-1)^m_., z_, n_, rm_, rp_, s_, op_, rest___ ] :>
  989.     Block [ {},
  990.         Assuming[ m < 0, dialogueAllFlag ];
  991.         If [ dialogueAllFlag, Print[ "where ", m, " is an integer" ] ];
  992.         Difference[m,n][ myinvz[x, z, n, rm, rp, s, op, rest] ] ] /;
  993.     FreeQ[m,n] && Implies[NumberQ[m], IntegerQ[m] && m > 0],
  994.  
  995.  
  996.   (*      B.  Summation                            *)
  997.   myinvz[ x_ / ( 1 - z_^-1 ), z_, n_, rest___ ] :>
  998.     Block [    {context, index},
  999.         context = $Context;
  1000.         $Context = "Global`";
  1001.         index = Unique["index"];
  1002.         $Context = context;
  1003.         Summation[index,0,n,1][ myinvz[x, z, index, rest] ] ],
  1004.  
  1005.  
  1006.   (*      C.  Periodic sequence with period k.                *)
  1007.   myinvz[ x_ / ( 1 - z_^k_), z_, n_, rm_, rp_, s_, op_, rest___ ] :>
  1008.     Block [ {},
  1009.         Assuming[ k < 0, dialogueAllFlag ];
  1010.         If [ dialogueAllFlag, Print[ "where ", k, " is an integer" ] ];
  1011.         Periodic[-k,n][ myinvz[x, z, n, rm, rp, s, op, rest] ] ]  /;
  1012.     FreeQ[k,n] && Implies[ NumberQ[k], IntegerQ[k] && ( k < 0 ) ],
  1013.  
  1014.  
  1015.  
  1016.  
  1017.   (*  V.  Z - T R A N S F O R M    S T R A T E G I E S            *)
  1018.  
  1019.   (*        If the following three rules are considered, then the    *)
  1020.   (*  the current expression could not be inverted by table lookup    *)
  1021.   (*  and applying properties.  So, the following strategy is used:    *)
  1022.  
  1023.   (*      A.  Factor terms inside of a log or exponential function.    *)
  1024.   (*      B.  Normalize the denominator.                *)
  1025.   (*      C.  Expand every term in the z-transform to its simplest    *)
  1026.   (*          components and try again.                    *)
  1027.   (*      D.  Substitute 1/z for z, take the inverse z-transform, and    *)
  1028.   (*          substitute -n for n in the result.            *)
  1029.   (*      E.  Complex cepstrum.                        *)
  1030.   (*      F.  If the current expression is a rational polynomial,    *)
  1031.   (*          then return an ARMA signal processing expression since    *)
  1032.   (*          it could not be inverted using partial fractions.        *)
  1033.   (*      G.  Perform a Laurent series expansion using the first    *)
  1034.   (*          Terms terms (Terms is an option) and try again.        *)
  1035.  
  1036.  
  1037.   (*      A.  Factor terms inside of a logarithmic or exponential expr.    *)
  1038.   myinvz[ Log[x_], z_, n_, rm_, rp_, s_, rest___ ] :>
  1039.     myinvz [ Log[Factor[x]], z, n,
  1040.          rm, rp, SetStateField[s, logformfield, False], rest ] /;
  1041.     GetStateField[s, logformfield] && PolynomialQ[x, z],
  1042.  
  1043.   myinvz[ x_. b_^(exp_), z_, n_, rm_, rp_, s_, rest___ ] :>
  1044.     myinvz [ x b^(Factor[exp]), z, n,
  1045.          rm, rp, SetStateField[s, expformfield, False], rest ] /;
  1046.     GetStateField[s, expformfield] && FreeQ[b, z],
  1047.  
  1048.  
  1049.   (*      B.  Expand definitions of all expressions to simplest form    *)
  1050.   myinvz[ x_, z_, n_, rm_, rp_, s_, rest___ ] :>
  1051.     myinvz [ ExpandAll[x], z, n, rm, rp,
  1052.          SetStateField[s, expandallfield, False], rest ] /;
  1053.     GetStateField[s, expandallfield],
  1054.  
  1055.  
  1056.   (*      C.  Substitute z^-1 for z, take the inverse z-transform,    *)
  1057.   (*          and substitute -n for n in the result.            *)
  1058.   myinvz[ x_, z_, n_, rm_, rp_, s_, rest___ ] :>
  1059.     Block [    {newrm, newrp, newx, seq},
  1060.         newx = x /. z -> z^-1;
  1061.         newrp = If [ ZeroQ[rm], Infinity, 1/rm ];
  1062.         newrm = If [ InfinityQ[rp], 0, 1/rp ];
  1063.  
  1064.         seq = MyInvZTransform [ newx, z, n, newrm, newrp,
  1065.                     anticausalInvZstate[], rest ];
  1066.         If [ InvalidInvZTransformQ[seq],
  1067.              myinvz [ x, z, n, rm, rp,
  1068.                   SetStateField[s, anticausalfield, False], rest ],
  1069.              seq /. n -> -n ] ] /;
  1070.     GetStateField[s, anticausalfield],
  1071.  
  1072.  
  1073.   (*      D.  Complex cepstrum [O&S, 498]                *)
  1074.   myinvz[ Log[x_], z_, n_, rm_, rp_, s_, rest___ ] :>
  1075.     -1/n myinvz [ z D[x, z] / x, z, n, rm, rp,
  1076.               SetStateField[s, cepstrumfield, False], rest ] /;
  1077.     GetStateField[s, cepstrumfield] &&
  1078.       FreeQ[ D[x, z], Derivative[1][f_][z] ],
  1079.  
  1080.  
  1081.   (*      E.  Invert rational polynomial as a cascade of an        *)
  1082.   (*          FIR filter and an IIR filter.                *)
  1083.   myinvz[ x_, z_, n_, rm_, rp_, s_, rest___ ] :>
  1084.     Block [    {firlist, iirlist, maxexp, zdenom, znumer},
  1085.         firlist = CoefficientList[znumer /. z -> z^-1, z];
  1086.         iirlist = CoefficientList[zdenom /. z -> z^-1, z];
  1087.         If [ SameQ[znumer, 1],
  1088.              Shift[ maxexp, n ]
  1089.               [ IIR[n, iirlist] ],
  1090.              Shift[ maxexp, n ]
  1091.               [ FIR[n, firlist]
  1092.                 [ IIR[n, iirlist] ] ] ] ] /;
  1093.     RationalFunctionQ[x, z] && RationalPolynomialQ[x, z]
  1094.  
  1095. }
  1096.  
  1097.  
  1098. (*  E N D     P A C K A G E  *)
  1099.  
  1100. End[]
  1101. EndPackage[]
  1102.  
  1103. If [ TrueQ[ $VersionNumber >= 2.0 ],
  1104.      On[ General::spell ];
  1105.      On[ General::spell1 ] ];
  1106.  
  1107.  
  1108. (*  H E L P     I N F O R M A T I O N  *)
  1109.  
  1110. Combine[ SPfunctions, {InvZTransform} ]
  1111. Protect[ InvZTransform ]
  1112.  
  1113.  
  1114. (*  E N D I N G     M E S S A G E  *)
  1115.  
  1116. Print["The inverse z-transform rules have been loaded."]
  1117. Null
  1118.